#Project: TV polarization
#Author: Josh Ryan and Jeff Lyons
#Working Date: Winter 23
#Task: Use Imai, Kim, Wang PanelMatch analysis--find best matching algorithm. Histogram of matched sets is in analysis file
#Related files: chamber_level_forstata_v4_102524.dta is loaded


library(foreign)
library(effects)
library(sciplot)
library(sandwich)
library(lmtest)   
library(BMS)
library(rJava)
library(lattice)
library(plyr)
library(lme4)
library(gmodels)
library(MASS)
library(stringr)
library("XLConnect")
library(reshape2)
library(xlsx)
library(dplyr)
library(doBy)
library(Hmisc)
library(stargazer)
library(xtable)
library(dplyr)
library(xtable)
library(foreign)
library(dotwhisker)#https://cran.r-project.org/web/packages/dotwhisker/vignettes/dotwhisker-vignette.html
library(berryFunctions) #for the insert row function
library(ggplot2)
library(gridExtra)
library(wnominate)
library(tidyr)
library(wnominate)
library(oc)
library(haven)
library(PanelMatch)

dat<-read_dta("C:/Users/A02177653/Dropbox/TV cameras-polarization/Data/Scripts/Replication/chamber_level_forstata_v4_102524.dta")

#Have to keep making as dataframe because of the way Panelmatch works
dat<-as.data.frame(dat)

#create state-chamber id that is character to use as labels--also because Panelmatch is picky
dat$state_cham2<-paste(dat$state, dat$Chamber,sep = "-")



#####Figures E1-E6 Covariate Balancing for Different Methods

####Late Budget
#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.latebudget<-dat[c("late_budget", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.latebudget<-as.data.frame(dat.latebudget)
dat.latebudget$state_cham<-as.integer(dat.latebudget$state_cham)
dat.latebudget$year<-as.integer(dat.latebudget$year)

#check balance improvement using late budget as outcome. first run matching results, then look at balance improvement
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "none", 
                         data = dat.latebudget, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                                  I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "late_budget",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.latebudget)
summary(PE.results)
plot(PE.results)

results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                           treatment = "treatment", refinement.method = "mahalanobis", 
                           data = dat.latebudget, match.missing = TRUE, 
                           covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "late_budget",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.latebudget)
summary(PE.results)
plot(PE.results)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                          treatment = "treatment", refinement.method = "ps.match", 
                          data = dat.latebudget, match.missing = TRUE, 
                          covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "late_budget",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.latebudget)
summary(PE.results)
plot(PE.results)

results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "ps.weight", 
                         data = dat.latebudget, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "late_budget",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.latebudget)
summary(PE.results)
plot(PE.results)

#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
par(mfrow=c(1,3))
balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.mah$att),
                data = dat.latebudget,
                main = "Mahalanobis Distance",
                covariates = c("late_budget", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.ps$att),
                data = dat.latebudget,
                main = "Propensity Score",
                covariates = c("late_budget", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.weight$att),
                data = dat.latebudget,
                main = "Propensity Score Weighting",
                covariates = c("late_budget", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

#confirmed by this test--pretty close between ps and weight, but ps wins

get_covariate_balance(results.none$att, dat.latebudget, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                          "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.mah$att, dat.latebudget, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                            "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.latebudget, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                          "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.latebudget, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                           "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

###############ps weight performs the best



###now do with kurtosis



#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.kurtosis<-dat[c("budget_kurtosis", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.kurtosis<-as.data.frame(dat.kurtosis)
dat.kurtosis$state_cham<-as.integer(dat.kurtosis$state_cham)
dat.kurtosis$year<-as.integer(dat.kurtosis$year)

#check balance improvement using late budget as outcome. first run matching results, then look at balance improvement
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.kurtosis, match.missing = TRUE, 
                           covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "budget_kurtosis",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.kurtosis)
summary(PE.results)
plot(PE.results)

results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.kurtosis, match.missing = TRUE, 
                          covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "budget_kurtosis",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.kurtosis)
summary(PE.results)
plot(PE.results)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.kurtosis, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "budget_kurtosis",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.kurtosis)
summary(PE.results)
plot(PE.results)

results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.kurtosis, match.missing = TRUE, 
                             covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "budget_kurtosis",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.kurtosis)
summary(PE.results)
plot(PE.results)

#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
par(mfrow=c(1,3))
balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.mah$att),
                data = dat.kurtosis,
                main = "Mahalanobis Distance",
                covariates = c("budget_kurtosis", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.ps$att),
                data = dat.kurtosis,
                main = "Propensity Score",
                covariates = c("budget_kurtosis", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.weight$att),
                data = dat.kurtosis,
                main = "Propensity Score Weighting",
                covariates = c("budget_kurtosis", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

#confirmed by this test--pretty close between ps and weight, but ps wins

get_covariate_balance(results.none$att, dat.kurtosis, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                       "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.mah$att, dat.kurtosis, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                      "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.kurtosis, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                     "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.kurtosis, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                         "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))


###############ps again performs the best





#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.diffs<-dat[c("diffs", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.diffs<-as.data.frame(dat.diffs)
dat.diffs$state_cham<-as.integer(dat.diffs$state_cham)
dat.diffs$year<-as.integer(dat.diffs$year)

#check balance improvement using late budget as outcome. first run matching results, then look at balance improvement
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.diffs, match.missing = TRUE, 
                           covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "diffs",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.diffs)
summary(PE.results)
plot(PE.results)

results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.diffs, match.missing = TRUE, 
                          covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "diffs",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.diffs)
summary(PE.results)
plot(PE.results)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.diffs, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "diffs",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.diffs)
summary(PE.results)
plot(PE.results)

results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.diffs, match.missing = TRUE, 
                             covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "diffs",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.diffs)
summary(PE.results)
plot(PE.results)

#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
par(mfrow=c(1,3))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.mah$att),
                data = dat.diffs,
                main = "Mahalanobis Distance",
                covariates = c("diffs", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.ps$att),
                data = dat.diffs,
                main = "Propensity Score",
                covariates = c("diffs", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.weight$att),
                data = dat.diffs,
                main = "Propensity Score Weighting",
                covariates = c("diffs", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

#confirmed by this test--pretty close between ps and weight, but ps wins

get_covariate_balance(results.none$att, dat.diffs, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                     "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.mah$att, dat.diffs, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                    "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.diffs, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                   "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.diffs, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                       "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))


###############ps again performs the best






#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.dem_sd<-dat[c("dem_sd", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.dem_sd<-as.data.frame(dat.dem_sd)
dat.dem_sd$state_cham<-as.integer(dat.dem_sd$state_cham)
dat.dem_sd$year<-as.integer(dat.dem_sd$year)

#check balance improvement using late budget as outcome. first run matching results, then look at balance improvement
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.dem_sd, match.missing = TRUE, 
                           covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "dem_sd",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.dem_sd)
summary(PE.results)
plot(PE.results)

results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.dem_sd, match.missing = TRUE, 
                          covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "dem_sd",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.dem_sd)
summary(PE.results)
plot(PE.results)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.dem_sd, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "dem_sd",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.dem_sd)
summary(PE.results)
plot(PE.results)

results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.dem_sd, match.missing = TRUE, 
                             covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "dem_sd",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.dem_sd)
summary(PE.results)
plot(PE.results)

#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
par(mfrow=c(1,3))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.mah$att),
                data = dat.dem_sd,
                main = "Mahalanobis Distance",
                covariates = c("dem_sd", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.ps$att),
                data = dat.dem_sd,
                main = "Propensity Score",
                covariates = c("dem_sd", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.weight$att),
                data = dat.dem_sd,
                main = "Propensity Score Weighting",
                covariates = c("dem_sd", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

#confirmed by this test--pretty close between ps and weight, but ps wins

get_covariate_balance(results.none$att, dat.dem_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                  "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.mah$att, dat.dem_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                 "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.dem_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.dem_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                    "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))



#########rep sd.




#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.rep_sd<-dat[c("rep_sd", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.rep_sd<-as.data.frame(dat.rep_sd)
dat.rep_sd$state_cham<-as.integer(dat.rep_sd$state_cham)
dat.rep_sd$year<-as.integer(dat.rep_sd$year)

#check balance improvement using late budget as outcome. first run matching results, then look at balance improvement
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.rep_sd, match.missing = TRUE, 
                           covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "rep_sd",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.rep_sd)
summary(PE.results)
plot(PE.results)

results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.rep_sd, match.missing = TRUE, 
                          covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "rep_sd",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.rep_sd)
summary(PE.results)
plot(PE.results)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.rep_sd, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "rep_sd",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.rep_sd)
summary(PE.results)
plot(PE.results)

results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.rep_sd, match.missing = TRUE, 
                             covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "rep_sd",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.rep_sd)
summary(PE.results)
plot(PE.results)

#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
par(mfrow=c(1,3))
balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.mah$att),
                data = dat.rep_sd,
                main = "Mahalanobis Distance",
                covariates = c("rep_sd", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.ps$att),
                data = dat.rep_sd,
                main = "Propensity Score",
                covariates = c("rep_sd", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.weight$att),
                data = dat.rep_sd,
                main = "Propensity Score Weighting",
                covariates = c("rep_sd", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

#confirmed by this test--pretty close between ps and weight, but ps wins

get_covariate_balance(results.none$att, dat.rep_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                   "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.mah$att, dat.rep_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                  "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.rep_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                 "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.rep_sd, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                     "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))




#########leg. prod



#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.billprop<-dat[c("billprop", "treatment", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "state_cham")]

#have to constantly be making these as integers, super annoying
dat.billprop<-as.data.frame(dat.billprop)
dat.billprop$state_cham<-as.integer(dat.billprop$state_cham)
dat.billprop$year<-as.integer(dat.billprop$year)

#check balance improvement using late budget as outcome. first run matching results, then look at balance improvement
results.none <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                           treatment = "treatment", refinement.method = "none", 
                           data = dat.billprop, match.missing = TRUE, 
                           covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                             I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                           + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                             I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                           size.match = 10, qoi = "att" ,outcome.var = "billprop",
                           lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.none, data = dat.billprop)
summary(PE.results)
plot(PE.results)

results.mah <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                          treatment = "treatment", refinement.method = "mahalanobis", 
                          data = dat.billprop, match.missing = TRUE, 
                          covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                            I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                          + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                            I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                          size.match = 10, qoi = "att" ,outcome.var = "billprop",
                          lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.mah, data = dat.billprop)
summary(PE.results)
plot(PE.results)

results.ps <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                         treatment = "treatment", refinement.method = "ps.match", 
                         data = dat.billprop, match.missing = TRUE, 
                         covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                           I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                         + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                           I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                         size.match = 10, qoi = "att" ,outcome.var = "billprop",
                         lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.billprop)
summary(PE.results)
plot(PE.results)

results.weight <- PanelMatch(lag = 4, time.id = "year", unit.id = "state_cham", 
                             treatment = "treatment", refinement.method = "ps.weight", 
                             data = dat.billprop, match.missing = TRUE, 
                             covs.formula = ~ I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
                               I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4))
                             + I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
                               I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
                             size.match = 10, qoi = "att" ,outcome.var = "billprop",
                             lead = 0:4, forbid.treatment.reversal = TRUE)

PE.results <- PanelEstimate(sets = results.ps, data = dat.billprop)
summary(PE.results)
plot(PE.results)

#A dot below the 45 degree line implies that the standardized mean balance is improved after the refinement for a particular time-varying covariate.
par(mfrow=c(1,3))
balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.mah$att),
                data = dat.billprop,
                main = "Mahalanobis Distance",
                covariates = c("billprop", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.ps$att),
                data = dat.billprop,
                main = "Propensity Score",
                covariates = c("billprop", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

balance_scatter(non_refined_set = results.none$att,
                matched_set_list = list(results.weight$att),
                data = dat.billprop,
                main = "Propensity Score Weighting",
                covariates = c("billprop", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                               "citizen_ideo", "lnpop", "total_intro"))

#confirmed by this test--pretty close between ps and weight, but ps wins

get_covariate_balance(results.none$att, dat.billprop, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                   "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.mah$att, dat.billprop, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                  "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.ps$att, dat.billprop, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                 "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))

get_covariate_balance(results.weight$att, dat.billprop, covariates = c("veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp",
                                                                     "citizen_ideo", "lnpop", "total_intro"),
                      ylim = c(-2,2))


